home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / FACTR.FOR < prev    next >
Text File  |  1985-11-29  |  3KB  |  128 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE FACTR
  5. C
  6. C        PURPOSE
  7. C           FACTORIZATION OF THE MATRIX A INTO A PRODUCT OF A LOWER
  8. C           TRIANGULAR MATRIX L AND AN UPPER TRIANGULAR MATRIX U.  L HAS
  9. C           UNIT DIAGONAL WHICH IS NOT STORED.
  10. C
  11. C        USAGE
  12. C           CALL FACTR(A,PER,N,IA,IER)
  13. C
  14. C        DESCRIPTION OF PARAMETERS
  15. C           A      MATRIX A
  16. C           PER    ONE DIMENSIONAL ARRAY WHERE PERMUTATIONS OF ROWS OF
  17. C                  THE MATRIX ARE STORED
  18. C                  DIMENSION OF PER MUST BE GREATER THAN OR EQUAL TO N
  19. C           N      ORDER OF THE MATRIX A
  20. C           IA     SIZE OF THE FIRST DIMENSION ASSIGNED TO THE ARRAY A
  21. C                  IN THE CALLING PROGRAM WHEN THE MATRIX IS IN DOUBLE
  22. C                  SUBSCRIPTED DATA STORAGE MODE.  IA=N WHEN THE MATRIX
  23. C                  IS IN SSP VECTOR STORAGE MODE.
  24. C           IER    ERROR INDICATOR WHICH IS ZERO IF THERE IS NO ERROR,
  25. C                  AND IS THREE IF THE PROCEDURE FAILS.
  26. C
  27. C        REMARKS
  28. C           THE ORIGINAL MATRIX, A,IS REPLACED BY THE TRIANGULAR FACTORS
  29. C
  30. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  31. C           NONE
  32. C
  33. C        METHOD
  34. C           SUCCESSIVE COMPUTATION OF THE COLUMNS OF L AND THE
  35. C           CORRESPONDING ROWS OF U.
  36. C
  37. C        REFERENCES
  38. C           J. H. WILKINSON - THE ALGEBRAIC EIGENVALUE PROBLEM -
  39. C           CLARENDON PRESS, OXFORD, 1965. H. J. BOWDLER, R. S. MARTIN,
  40. C           G. PETERS, AND J. H. WILKINSON - 'SOLUTION OF REAL AND
  41. C           COMPLEX SYSTEMS OF LINEAR EQUATIONS', NUMERISCHE MATHEMATIK,
  42. C           VOL. 8, NO. 3, 1966, P. 217-234.
  43. C
  44. C     ..................................................................
  45. C
  46.       SUBROUTINE FACTR(A,PER,N,IA,IER)
  47.       DIMENSION A(1),PER(1)
  48.       DOUBLE PRECISION DP
  49. C
  50. C        COMPUTATION OF WEIGHTS FOR EQUILIBRATION
  51. C
  52.       DO 20 I=1,N
  53.       X=0.
  54.       IJ=I
  55.       DO 10 J=1,N
  56.       IF (ABS(A(IJ))-X)10,10,5
  57.     5 X=ABS(A(IJ))
  58.    10 IJ=IJ+IA
  59.       IF (X) 110,110,20
  60.    20 PER(I)=1./X
  61.       I0=0
  62.       DO 100 I=1,N
  63.       IM1=I-1
  64.       IP1=I+1
  65.       IPIVOT=I
  66.       X=0.
  67. C
  68. C        COMPUTATION OF THE ITH COLUMN OF L
  69. C
  70.       DO 50 K=I,N
  71.       KI=I0+K
  72.       DP=A(KI)
  73.       IF (I-1) 110,40,25
  74.    25 KJ=K
  75.       DO 30 J=1,IM1
  76.       IJ=I0+J
  77.       DP=DP-1.D0*A(KJ)*A(IJ)
  78.    30 KJ=KJ+IA
  79.       A(KI)=DP
  80. C
  81. C        SEARCH FOR EQUILIBRATED PIVOT
  82. C
  83.    40 IF (X-DABS(DP)*PER(K))45,50,50
  84.    45 IPIVOT=K
  85.       X=DABS(DP)*PER(K)
  86.    50 CONTINUE
  87.       IF (X)110,110,55
  88. C
  89. C        PERMUTATION OF ROWS IF REQUIRED
  90. C
  91.    55 IF (IPIVOT-I) 110,70,57
  92.    57 KI=IPIVOT
  93.       IJ=I
  94.       DO 60 J=1,N
  95.       X=A(IJ)
  96.       A(IJ)=A(KI)
  97.       A(KI)=X
  98.       KI=KI+IA
  99.    60 IJ=IJ+IA
  100.       PER(IPIVOT)=PER(I)
  101.    70 PER(I)=IPIVOT
  102.       IF (I-N) 72,100,100
  103.    72 IJ=I0+I
  104.       X=A(IJ)
  105. C
  106. C        COMPUTATION OF THE ITH ROW OF U
  107. C
  108.       K0=I0+IA
  109.       DO 90 K=IP1,N
  110.       KI=I0+K
  111.       A(KI)=A(KI)/X
  112.       IF (I-1)110,90,75
  113.    75 IJ=I
  114.       KI=K0+I
  115.       DP=A(KI)
  116.       DO 80 J=1,IM1
  117.       KJ=K0+J
  118.       DP=DP-1.D0*A(IJ)*A(KJ)
  119.    80 IJ=IJ+IA
  120.       A(KI)=DP
  121.    90 K0=K0+IA
  122.   100 I0=I0+IA
  123.       IER=0
  124.       RETURN
  125.   110 IER=3
  126.       RETURN
  127.       END
  128.